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Abstract. We consider the problem of characterisation of burst sources detected with the 
Laser Interferometer Space Antenna (LISA) using the multi-modal nested sampling algorithm, 
MULTlNEST. We use MULTlNEST as a tool to search for modelled bursts from cosmic 
string cusps, and compute the Bayesian evidence associated with the cosmic string model. 
As an alternative burst model, we consider sine-Gaussian burst signals, and show how the 
evidence ratio can be used to choose between these two alternatives. We present results from 
an application of MULTlNEST to the last round of the Mock LISA Data Challenge, in which 
we were able to successfully detect and characterise all three of the cosmic string burst sources 
present in the release data set. We also present results of independent trials and show that 
MULTlNEST can detect cosmic string signals with signal-to-noise ratio (SNR) as low as ~ 7 
and sine-Gaussian signals with SNR as low as ~ 8. In both cases, we show that the threshold 
at which the sources become detectable coincides with the SNR at which the evidence ratio 
begins to favour the correct model over the alternative. 
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1. Introduction 

The development of data analysis algorithms for gravitational wave detectors on the ground 
and for the proposed space-based gravitational wave detector, the Laser Interferometer Space 
Antenna (LISA) [ 1 1, is an active area of current research. LISA data analysis activities are 
being encouraged by the Mock LISA Data Challenges [2 1 (MLDCs). Prior to the most recent 
round, round 3, which finished in April 2009, the MLDCs had included data sets containing 
individual and multiple white-dwarf binaries, a realisation of the whole galaxy of compact 
binaries, single and multiple non-spinning supermassive black hole (SMBH) binaries, either 
isolated or on top of a galactic confusion background and isolated extreme-mass-ratio 
inspiral (EMRJ) sources in purely instrumental noise. Round 3 included a realisation of 
the galactic binaries that included chirping systems, a single data set containing multiple 
overlapping EMRIs and data sets containing three new types of source — inspirals of spinning 
supermassive black hole binaries, bursts from cosmic string cusps and a stochastic background 
of gravitational radiation [2|. Several of these sources, including the EMRI signals, the 
spinning black hole binaries and the cosmic string bursts, have highly multi-modal likelihood 
surfaces which can cause problems for algorithms such as Markov Chain Monte Carlo 
(MCMC), as these may become stuck in secondary maxima rather than finding the primary 
mode J5] 131 15] |6) . Recently [7| we applied a new algorithm to the problem of gravitational 
wave data analysis, MultiNest 18,91, which is a nested sampling algorithm, optimized for 
problems with highly multi-modal posteriors. We demonstrated that MultiNest could be 
used as a search tool and to recover posterior probability distributions for the case of data 
sets containing multiple signals from non-spinning SMBH binary inspirals. We used the 
algorithm as a search tool in MLDC round 3 to analyse two of the data sets — the spinning 
SMBH inspirals (challenge 3.2) and the cosmic string bursts (challenge 3.4). We will discuss 
the second of these applications in this paper. 

The nested sampling algorithm [ 1 1 was developed as a tool for evaluating the Bayesian 
evidence. It employs a set of live points, each of which represents a particular set of 
parameters in the multi-dimensional search space. These points move as the algorithm 
progresses, and climb together through nested contours of increasing likelihood. At each step 
the algorithm works by finding a point of higher likelihood than the lowest likelihood point 
in the live point set and then replacing the lowest likelihood point with the new point. The 
difficulty is to sample, efficiently, new points of higher likelihood from the remaining prior 
volume. MultiNest achieves this using an ellipsoidal rejection sampling scheme. It has 
demonstrated orders of magnitude improvement over standard methods in several applications 
in cosmology and particle physics ifTTl [121 [131 [14*1 . It also proved to be very effective in a 
gravitational wave context for the test case mentioned above Q. 

MultiNest can be used as a search tool and returns the posterior probability 
distributions as a by-product. We employed the algorithm successfully to analyse MLDC 
challenge 3.4, accurately finding all three of the cosmic string burst sources present in the 



data set. This will be discussed in more detail in Section 3.1 The cosmic string cusp is a 
special type of burst source, since the waveform is known and can be modelled, so a matched 
filtering search is possible and that is what MultiNest does in computing the likelihoods 
across the parameter space. However, there may be bursts of gravitational waves in the LISA 
data stream from other sources and the question arises as to whether we would be able to 
detect them and if we can distinguish cosmic string bursts from bursts due to these other 
sources. The evidence that MultiNest computes is one tool that can be used to address 
model selection. Evidence has been used for gravitational wave model selection to test the 
theory of relativity in ground based observations [ 15 1, and, in a LISA context, to distinguish 
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between an empty data set and one containing a signal |[T6l . Calculation of an evidence ratio 
requires a second model for the burst as an alternative to compare against. We chose to use a 
sine-Gaussian waveform as a generic alternative burst model, as this is one of the burst models 
commonly used in LIGO data analysis (see for instance ifTTl ). We find that the evidence is a 
powerful tool for characterising bursts — for the majority of detectable bursts, the evidence 
ratio strongly favours the true model over the alternative. While the sine-Gaussian model 
does not necessarily well describe all possible un-modelled bursts, these results suggest that 
the evidence can be used to correctly identify any cosmic string bursts that are present in the 
LISA data stream. 

This paper is organised as follows. In Section [2] we describe the methods that we 
employ in this analysis — we describe Bayesian inference, nested sampling, the MultiNest 
algorithm, the two burst waveform models we use and the detector noise spectral density. 
In Section [5] we present our results, including a description of our methods and results for 
the analysis of MLDC challenge 3.4 (in Section 3.1 1; an estimate of the signal-to-noise ratio 
(SNR) required for detection of cosmic string and sine-Gaussian bursts using the MultiNest 
algorithm; and a calculation of the evidence ratio of the two models for data sets containing 
each type of source. We finish in Section |4] with a summary and discussion of directions for 
further research. 



2. Methodology 

2.1. Bayesian Inference 

Bayesian inference methods provide a consistent approach to the estimation of a set of 
parameters in a model (or hypothesis) H for the data D. Bayes' theorem states that 

where Pr(0|D, H) = P{&) is the posterior probability distribution of the parameters, 
Pr(D|0,if) = £(0) is the likelihood, Pv(&\H) = tt(0) is the prior distribution, and 
Pr(D|if) = Z is the Bayesian evidence. 

Bayesian evidence is the factor required to normalise the posterior over 0: 

Z = J £(®)n{@)d N ®, (2) 

where N is the dimensionality of the parameter space. Since the Bayesian evidence is 
independent of the parameter values 0, it is usually ignored in parameter estimation problems 
and the posterior inferences are obtained by exploring the un-normalized posterior using 
standard MCMC sampling methods. 

Bayesian parameter estimation has been used quite extensively in a variety of 
astronomical applications, including gravitational wave astronomy, although standard MCMC 
methods, such as the basic Metropolis-Hastings algorithm or the Hamiltonian sampling 
technique (see e.g. 1 18 1), can experience problems in sampling efficiently from a multi-modal 
posterior distribution or one with large (curving) degeneracies between parameters. Moreover, 
MCMC methods often require careful tuning of the proposal distribution to sample efficiently, 
and testing for convergence can be problematic. 

In order to select between two models Hq and Hi one needs to compare their respective 
posterior probabilities given the observed data set D, as follows: 

Pr(i2i|D) Pr(D|i?i)Pr(iJi) Z 1 Pr{H 1 ) 



Pr(tf |D) Pr(D|ff )Pr(iJ ) Z Py(H ) ' 



(3) 
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where Pr(iJi)/ Vi(Hq) is the prior probability ratio for the two models, which can often 
be set to unity but occasionally requires further consideration (see, for example, ifTTl fl2l for 
some examples where the prior probability ratio should not be set to unity). It can be seen from 
Eq. (|3]l that the Bayesian evidence plays a central role in Bayesian model selection. As the 
average of likelihood over the prior, the evidence automatically implements Occam's razor: 
a simpler theory which agrees well enough with the empirical evidence is preferred. A more 
complicated theory will only have a higher evidence if it is significantly better at explaining 
the data than a simpler theory. 

Unfortunately, evaluation of Bayesian evidence involves the multidimensional integral 
(Eq. (|2]l) and thus presents a challenging numerical task. Standard techniques like 
thermodynamic integration lfT9l are extremely computationally expensive which makes 
evidence evaluation typically at least an order of magnitude more costly than parameter 
estimation. Some fast approximate methods have been used for evidence evaluation, such 
as treating the posterior as a multivariate Gaussian centred at its peak (see, for example, 
[20 1), but this approximation is clearly a poor one for highly non-Gaussian and multi-modal 
posteriors. Various alternative information criteria for model selection are discussed in ll2D . 
but the evidence remains the preferred method. 



2.2. Nested Sampling and the MULTlNEST Algorithm 

Nested sampling iflOl is a Monte Carlo method targetted at the efficient calculation of the 
evidence, but also produces posterior inferences as a by-product. It calculates the evidence 
by transforming the multi-dimensional evidence integral into a one-dimensional integral that 
is easy to evaluate numerically. This is accomplished by defining the prior volume X as 
dX = Tr(@)d D @, so that 

X(X)= f n(®)d N ®, (4) 

JC(&)>\ 

where the integral extends over the region(s) of parameter space contained within the iso- 
likelihood contour £(@) = A. The evidence integral, Eq. (j2j), can then be written as 

Z = [ £{X)dX, (5) 







where C(X), the inverse of Eq. is a monotonically decreasing function of X. Thus, if one 
can evaluate the likelihoods Ci = £(Xi), where Xj is a sequence of decreasing values, 

< X M < ■ ■ ■ < X 2 < Xi < X = 1, (6) 

as shown schematically in Fig. [T] the evidence can be approximated numerically using 
standard quadrature methods as a weighted sum 

AI 

Z = ^dwi, (7) 

i=l 

where the weights Wi for the simple trapezium rule are given by Wi = ^(X^_i — Xi + i). An 
example of a posterior in two dimensions and its associated function C(X) is shown in Fig. |T| 

2.2.7. Evidence Evaluation The summation in Eq. (j7]l is performed as follows. The iteration 
counter is first set to i = and N 'active' (or 'live') samples are drawn from the full 
prior 7r(@), so the initial prior volume is Xq = 1. The samples are then sorted in order 
of their likelihood and the smallest (with likelihood Co) is removed from the active set (hence 
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becoming 'inactive') and replaced by a point drawn from the prior subject to the constraint 
that the point has a likelihood C > Cq. The corresponding prior volume contained within the 
iso-likelihood contour associated with the new live point will be a random variable given 
by X\ = t\X§, where t\ follows the distribution Pr(<) = Nt N ^ x (i.e., the probability 
distribution for the largest of N samples drawn uniformly from the interval [0, 1]). At each 
subsequent iteration i, the removal of the lowest likelihood point £, in the active set, the 
drawing of a replacement with C > Ci and the reduction of the corresponding prior volume 
X.- L = tiXi-i are repeated, until the entire prior volume has been traversed. The algorithm 
thus travels through nested shells of likelihood as the prior volume is reduced. The mean and 
standard deviation of log f , which dominates the geometrical exploration, are: 

E[\ogt] = -l/N, a[\ogt] = 1/N. (8) 

Since each value of log t is independent, after i iterations the prior volume will shrink down 
such that \ogXi rj —{i ± \/i)/N. Thus, one takes Xi — exp(—i/N). 



2.2.2. Stopping Criterion The nested sampling algorithm is terminated when the evidence 
has been computed to a pre-specified precision. The evidence that could be contributed by 
the remaining live points is estimated as AZj = £ max Xi, where £ max is the maximum- 
likelihood value of the remaining live points, and Xi is the remaining prior volume. The 
algorithm terminates when AZ; is less than a user-defined value (we use 0.5 in log-evidence). 



2.2.3. Posterior Inferences Once the evidence Z is found, posterior inferences can be easily 
generated using the final live points and the full sequence of discarded points from the nested 
sampling process, i.e., the points with the lowest likelihood value at each iteration i of the 
algorithm. Each such point is simply assigned the probability weight 

Pi = ^- (9) 

These samples can then be used to calculate inferences of posterior parameters such as 
means, standard deviations, covariances and so on, or to construct marginalised posterior 
distributions. 
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2.2.4. MultiNest Algorithm The most challenging task in implementing nested sampling 
is to draw samples from the prior within the hard constraint C > £, at each iteration i. The 
MultiNest algorithm [8,9] tackles this problem through an ellipsoidal rejection sampling 
scheme. The live point set is enclosed within a set of (possibly overlapping) ellipsoids and 
a new point is then drawn uniformly from the region enclosed by these ellipsoids. The 
ellipsoidal decomposition of the live point set is chosen to minimize the sum of volumes 
of the ellipsoids. The ellipsoidal decomposition is well suited to dealing with posteriors 
that have curving degeneracies, and allows mode identification in multi-modal posteriors. 
If there are subsets of the ellipsoid set that do not overlap with the remaining ellipsoids, these 
are identified as a distinct mode and subsequently evolved independently. The MultiNest 
algorithm has proven very useful for tackling inference problems in cosmology and particle 
physics |[TT1fT^[T3lfT4l . typically showing two orders of magnitude improvement in efficiency 
over conventional techniques. More recently, it has been shown to perform well as a search 
tool for gravitational wave data analysis Q. 



2.3. Burst waveform models 

2.3.1. Cosmic Strings Gravitational waves can be generated by cosmic strings through 
the formation of cusps, where portions of the string are traveling at nearly the speed of 
light [22, 23 1. Such radiation is highly beamed and, when viewed along the emission axis, 
is linearly polarized and takes a simple power-law form, hit) oc \t — t^ 1 / 3, 11221 . When 
viewed slightly off-axis, the waveform is still approximately linearly polarized but the cusp 
spectrum is rounded off and decays rapidly for frequencies above / max ~ 2/(a 3 L), where a 
is the viewing angle and L is the dimension of the feature generating the cusp [22, 24J. The 
particular model for the frequency domain waveform adopted in the MLDC is given by [2| 

, , , f <Af ^ I ^ /max 

1^)1 = 1^-^(1-^) />/_ ' ™ 

In addition, the MLDC waveforms include a fourth-order Butterworth filter to mitigate 
dynamic-range issues associated with inverse Fourier transforms. We adopt the same ansatz 
as the MLDC, namely that the Fourier domain waveform amplitude is given by 

-4 



\h+\ = 




f < fn 



-4 

exp f 1 - -J— ) / > f n 



\h x \=0 (11) 
and the phase by exp(27ri/t c ), where t c is the burst time. 



2.3.2. Sine Gaussians A sine-Gaussian waveform is centred on a particular frequency, and 
has exponentially suppressed power at nearby frequencies. We choose to consider a lineraly 
polarised sine-Gaussian, for which the waveform magnitudes in the frequency domain are 
given by 

M = 5iS^(t(t)} l "*' = ° (12) 

where A is the dimensionless amplitude, f c is the frequency of the sinusoidal oscillation, and 
Q is the width. The phase of the wave is again exp(27n/£ c ), where t c is the burst time. In the 
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time-domain, the sine-Gaussian is a small burst "packet" of particular frequency, f c , with the 
number of cycles in the burst determined by Q. 

2.4. Detector model 

To include the LISA response we made use of the static LISA model as described in |25|. 
This approximation is valid for burst sources, as LISA does not move significantly over the 
typically short duration of the bursts, which is of the order of 1000s. The static LISA response 
model was also adopted for cosmic string bursts in |24|. 

Three optimal detection time-delay interferometry (TDI) channels, A, E and T, can be 
constructed from the LISA data stream [26|. The noise is uncorrected within and across these 
channels. For the search of the MLDC data, the A, E and T channels were constructed as 
linear combinations of the three Michelson channels, X, Y and Z, that were provided in the 
data release. The power spectral densities of the noise in the three channels are given by 

S A (f) = S E (f) = 16 sin 2 (2^/i L ) (2 (l + cos(27r ft L ) + cos 2 (2^ L )) S pjn (f) 

+ (l + cos(2Trft L )/2)S sr J 2 ) (13) 

S T (f) = 16 sin 2 (2^ L ) (2(1-2 cos(2tt ft L ) + cos 2 (27r/t L )) S pm {f) 

+ (l-cos(2^/t L ))5 sn / 2 ) (14) 

v,</) =(i + (^) 2 )^ 

(15) 

where tz = 16.678s is the light travel time along one arm of the LISA constellation, 
•S'acc = 2.5 x 10 _48 Hz~ 1 is the proof mass acceleration noise and S sn = 1.8 x 10~ 37 Hz -1 
is the shot noise. 

The LISA data stream will also contain a confusion noise foreground from unresolved 
white dwarf binaries in our galaxy. These were not included in the noise model for the MLDC 
round 3.4, and so we have ignored the confusion noise contribution in the current work. In 
practice, we found that the T channel was too noisy to be used in the MLDC search, and so 
we used the A and E channels only for data analysis. 

2.5. Likelihood evaluation 

The space of gravitational waveform signals possesses a natural scalar product Il27ll28l 

(/i|s)=2 iTsS) [ jiif)r{f)+k{mf) }' (i6) 

where 

h(f) = / dth(t)e 27rlft (17) 



is the Fourier transform of the time domain waveform h(t). The quantity S n (f) is the one- 
sided noise spectral density of the detector. For a family of sources with waveforms h(t; A) 
that depend on parameters A, the output of the detector, s(t) = h(t; A ) + n(t), consists of 
the true signal h(t; A ) and a particular realisation of the noise, n(t). Assuming that the noise 
is stationary and Gaussian, the logarithm of the likelihood that the parameter values are given 
by A is 



log £ ( A ) =C-(s-h A 



MA)) /2, (18) 
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and it is this log-likelihood that is evaluated at each point in the search. The constant, C, 
depends on the dimensionality of the search space, but its value is not important as we are only 
interested in the relative likelihoods of different points. As mentioned earlier, the LISA data 
stream has several independent data channels. The total multi-channel likelihood is obtained 
by summing the scalar product ( 16 1 over the various channels (in our work we use only the 
two channels, A and E), 



2.6. Priors 

The parameter space over which to search for signals must also be specified. For the searches 
with the cosmic string model, we used uniform priors for each parameter that covered the 
signal space from which the MLDC sources were drawn. These correspond to the following 
ranges: 

log(A) G [-23.3, -21], / max G [0.001,0.5] Hz, 

t c e [T ,T + AT]s, sinOS) G [-1,1], 0€[O,27r], ^e[0,2ir]. 
Here, j3 is the sky latitude, <f> is the sky azimuthal angle, and i/) is the polarization. The data 
stream was divided into segments for the search and To is the start time of the particular data 
segment being searched and AT is that segment's length. 

For the sine-Gaussian model, we maintained the approach of using uniform prior ranges 
while making sure that all signals would fall within our specified bounds. The priors used 
were 

log(A) G [-22.3, -18], f c £ [0.001, 0.5] Hz, Q G [1, 10] 

« c e[T 0l T + AT] S , sin(/3) G [-1,1], ^[0,24 V € [0, 2tt] . 

T , AT, f3, <fi, and tp are the same physical quantities as in the cosmic string model; the prior 
on log (A) was modified in order to cover the same range of SNRs for the signals. 



3. Results 



3.1. MLDC challenge 3.4 search 

Mock LISA Data Challenge 3.4 (2] consisted of a data set of 2 21 samples at a cadence of 
Is, containing GW burst signals from cosmic string cusps [22, 23]. The signals occurred 
as a Poissonian random process throughout the data set, with an expectation value of five 
events. To analyse the data, we divided the data stream into segments of 32768s in length, 
with 2048s of overlap between neighbouring segments. The overlap was chosen to ensure that 
each of the bursts, which were of typical duration of ~ 1000s, would be entirely contained 
in at least one segment. To avoid artefacts in the Fourier domain, we used Welch windowing 
when performing the Fast Fourier Transform (FFT) of each segment. We also analysed data 
segments of the same length, offset from the first set by 16384s, to check no signals had been 
missed and to verify that the detected signals were not edge artefacts. We also repeated the 
search with segments of 16384s and 8192s in length to verify our results. MultiNest was 
run with 4000 live points on each segment, using 8 3GHz Intel Woodcrest processors in 
parallel; the algorithm took approximately 5 minutes to run on a single segment when no 
signal was present, and twice this when one was present. The search of the offset and shorter 
segments returned all of the same detections. We ran the search on the MLDC training data 
set to test the algorithm, and then on the challenge data set. 

To assess the quality of our solutions we will not only use the recovered parameters for a 
source, but also the signal-to-noise ratio (SNR) of the recovered solution. The degeneracies in 
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Figure 2. True signal and the nine modes of highest evidence identified by the algorithm in 
our search for the first source in the training data. 

the parameter space mean that a waveform can match the true waveform very well with very 
different parameter values. The recovered SNR is a measure of how much of a residual would 
be left in the data stream if the recovered parameters were used to subtract a given source from 
the data set. If the recovered SNR is close to the true SNR we can say that we have correctly 
identified the waveform of the source even if the parameters are quite different. 

3.1.1. Training data The MLDC training data set for challenge 3.4 contained five cosmic 
string bursts. Our search successfully identified the five sources, and recovered posterior 
probability distributions for their parameters. The posteriors were highly degenerate, 
containing many modes. These modes corresponded to degeneracies in the waveform 
parameter space that arise because the burst signals are of short duration compared to the 
orbital period of the detector. Such degeneracies were noted in |24|. We can use the 
local Bayesian evidence of the modes, as computed by MultiNest, to classify the various 
solutions. In Figure [2] we show the waveforms corresponding to the nine modes of highest 
local evidence identified in the search for the first of the signals in the training data set. We 
also show the true waveform signal for comparison. It is clear that all of the modes that were 
identified have nearly identical waveforms, so we would not expect to determine uniquely the 
true source parameters. However, it is also clear that MultiNest is successfully finding and 
identifying the signal. 

In Table [T] we summarise the parameters of the five signals in the training data set and 
the solutions recovered by MultiNest. For simplicity we present only the solution of 
highest local evidence identified by the algorithm, although in all cases we found multiple 
other modes. We quote both the maximum likelihood parameters found by MultiNest for 
this mode, and the mean and standard deviation computed from the posterior of that mode. 
We see that when the maximum frequency, f max , is significantly below Nyquist, then we can 
determine it well, but not when it is close to or above that frequency. This is entirely expected, 
as this parameter only modifies the signal for / > / max and therefore does not leave a 
signature in the data set when / max is too large. There is a degeneracy in the polarisation angle 
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corresponding to the shift ip ~ > tp + ir, and accounting for this the polarisation is moderately 
well determined, as is the amplitude. The time of coalesence at the detector is extremely well 
determined, but t c is the time of coalesence at the Solar System barycentre, and this is not 
so well determined due to uncertainties in the sky position. The three parameters, (3s,4>s, t c , 
show correlated degeneracies. These arise because the ability of LISA to resolve sky position 
relies on the motion of the detector around the Sun, but over the duration of a typical burst 
LISA is essentially static. At higher gravitational wave frequencies some of this degeneracy is 
weakly broken, but it is very severe at low frequency. There is also a well known degeneracy 
in the LISA response between antipodal sky positions, corresponding to gravitational waves 
crossing the detector in opposite directions. 

The standard deviations of the modes are quite large in all cases, which reflects the 
existence of these degeneracies. In most cases, the true parameters are within 3 standard 
deviations of the means, so the mean and error returned by MultiNest adequately reflect the 
parameters of the underlying signal. One interesting result is that, for source 4, the distribution 
of /max is as broad as it is in the three cases where this parameter could not be measured 
(sources 2, 3 and 5), and yet the maximum likelihood value of 0. 0118Hz is very close to the 
true value of 0. 0115Hz. For this source, the posterior recovered by MultiNest was very 
broad, with a narrow peak in the vicinity of the true value. This suggests that / max may be 
difficult to constrain even when it has a value significantly below Nyquist. 

The degeneracies in the parameter space are also evident in the marginalised posterior 
probability distributions, which we can construct from the MultiNest live point set 
evolution. In Figure [3] we show 2D and ID marginalised posterior distributions for the third 
of the sources present in the training data set. We note that the break frequency, / max , is 
unconstrained as it was only slightly below (0.464Hz) the Nyquist frequency of the data set 
(0.5Hz) in this case. The degeneracies in the posteriors for sky position and time of coalesence 
are also evident. We see two dominant peaks in the sky position posterior, which are points 
antipodal to one another on the sky. 

The existence of these degeneracies in the parameter space suggest that a comparison of 
recovered and true parameters is not the best way to characterise the quality of a detection, 
as mentioned above. An alternative approach is to compute the signal-to-noise ratio of the 
recovered source and compare that to the SNR of the injected source. This was possible 
using evaluation tools available in the LISATOOLS code librarjyj] In Table [2] we show the 
SNRs in the A and E channels for each of the five training sources. We give results for the 
recovered mode of highest local evidence, plus the corresponding result using the true source 
parameters. We see that, in all cases, the solution found by MultiNest recovered almost all 
of the SNR of the true source. If the burst source was subtracted from the data set using the 
recovered parameters, the residual in the data set would be at the level of instrumental noise 
fluctuations. We can therefore say that we correctly identified all of the sources. 

3.1.2. Challenge data Three signals were found in the challenge data set at t c ~ 6 x 10 5 s, 
t c ~ 1.07 x 10 6 s and t c ~ 1.6 x 10 6 s and this was the correct number of signals. We label 
these as source candidates 3, 2 and 1 respectively for consistency with the MLDC key file. For 
source candidates 1 and 3, MultiNest returned a flat distribution in / max , thereby leading 
us to conclude that the actual maximum was above the Nyquist frequency of 0.5 Hz. We 
identified several modes in the posterior in each case, and used the evidence to characterise 
these as for the training data. We decided to submit the two modes of highest evidence, as in 
the analysis of the training data these generally corresponded to two antipodal sky solutions. 

1 1 http://code.google.eom/p/lisatools/ 
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-0.553 


4.88 


0.24 


1.8892 x 10 6 




Std Dev 


±0.89 x 10~ 22 


±0.34 


±0.16 


±0.34 


±0.14 


±0.00014 x 10 6 




True 


6.60 x 10~ 22 


3.98 


-0.934 


0.677 


0.464 


1.8645 x 10 y 




Recovered 


9.96 x 10~ 22 


1.35 


-0.703 


1.55 


0.5* 


1.8644 x 10 6 


3 


Mean 


8.77 x 10~ 21 


1.34 


-0.835 


1.42 


0.26 


1.8644 x 10 6 




Std. Dev. 


±0.83 x 10" 22 


±0.19 


±0.14 


±0.22 


±0.14 


±0.000071 x 10 6 




True 


2.65 x 10~ 21 


3.54 


0.239 


1.09 


0.0115 


2.0603 x 10 6 




Recovered 


2.42 x 10~ 21 


0.505 


0.817 


2.35 


0.0118 


2.0598 x 10 6 


4 


Mean 


2.31 x 10~ 21 


0.517 


0.505 


2.50 


0.244 


2.0597 x 10 6 




Std. Dev. 


±0.30 x 10~ 21 


±1.3 


±0.20 


±0.10 


±0.15 


±0.000049 x 10 6 




True 


3.42 x 10~ 22 


1.09 


-1.03 


1.16 


2.28 


1.4983 x 10 6 




Recovered 


5.39 x 10~ 22 


4.50 


-0.838 


1.68 


0.5* 


1.4982 x 10 6 


5 


Mean 


4.23 x 10~ 22 


4.04 


-0.419 


2.99 


0.257 


1.4983 x 10 6 




Std. Dev. 


±1.1 x 10~ 22 


±1.0 


±0.51 


±2.01 


±0.14 


±0.00024 x 10 6 



Table 1. Parameters of the five signals present in the training data set, and the parameters 
recovered by MULTlNEST for the mode of highest local evidence identified by the search. 
The parameters are amplitude, A, polarisation, sky latitude, and longitude, 0g, break 
frequency, / max , and time of coalesence at the Solar System barycentre, t c . In the row labelled 
"Recovered" we quote the maximum likelihood parameters of the mode with highest evidence 
identified by MULTlNEST. In the rows labelled "Mean" and "Std. Dev." we quote the mean 
and standard deviation in that parameter, as computed from the recovered posterior for the 
mode of highest evidence. Where the value of /max is marked by a * , the distribution in this 
parameter was very flat, indicating that / max was probably above Nyquist. In these cases, we 
set /max to the Nyquist frequency of 0.5Hz for the maximum likelihood parameters, but still 
quote the actual recovered mean and standard deviation. 



Channel 


1 


2 


Source 
3 


4 


5 


A 


True 
Recovered 


52.4 
52.4 


12.9 
12.7 


16.1 
16.5 


12.6 
12.6 


4.34 
4.19 


E 


True 
Recovered 


17.1 
17.1 


17.7 
17.6 


26.8 
26.7 


76.6 
76.6 


13.6 
13.6 



Table 2. SNRs in the A and E channels for the five sources present in the training data, as 
computed using the evaluation scripts in LISATOOLS. We quote results using the true source 
parameters, labelled "True", and using the parameters of the mode identified by MULTlNEST 
of highest local evidence, labelled "Recovered". 
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Figure 3. Two-dimensional marginalised posteriors as recovered by MULTlNEST in the search 
for the third of the cosmic string bursts in the MLDC training data set. The parameters, from 
top-to-bottom and left-to-right, are colatitude, longitude, burst time, burst amplitude, burst 
break frequency and waveform polarization. At the top of each column we also show the 
one-dimensional posterior for the column parameter. 



However, for source candidate 3 there was a third mode of almost equal local evidence, and so 
we submitted that mode as well. The middle source, 2, was found to have a break frequency 
of 0.0011 Hz, very close to the minimum value in the prior range of 0.001 Hz. At this low 
frequency, LISA is not able to resolve the sky position for a burst source. As expected, we 
found very broad posteriors for all parameters other than the maximum frequency. As the 
posterior was not well separated into modes, we chose to submit only one set of parameters 
in this case, which we took to be the maximum a-posteriori parameters. 

In Table [3] we list the true parameters for the three sources in the challenge data sej^j] 
and the parameters recovered by MultiNest. For simplicity we only include the parameters 
for the mode of highest local evidence, and not all of the modes submitted. The agreement 

1 Available at http : / /www .tapir.caltech. edu/~mldc/ result s3 / 
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Source 


A 


i 


0s 


(j>S 


/may (HZ) 


t c (s) 




True 


8.66 x 10~ 22 


3.32 


0.556 


3.71 


0.0296 


1.6022 x 10 B 


1 


Recovered 


1.65 x lCT 21 


2.82 


0.349 


6.01 


0.5* 


1.6030 x 10 6 


Mean 


1.13 x lCT 21 


2.70 


-0.141 


4.38 


0.27 


1.6030 x 10 6 




Std. Dev. 


±0.22 x 10" 21 


±0.14 


±0.30 


±2.6 


±0.14 


±0.000071 x 10 6 




True 


2.79 x 10" 21 


5.12 


-0.444 


3.17 


0.00108 


1.0727 x 10 6 




Recovered 


2.97 x 10~ 21 


0.271 


-0.893 


5.12 


0.00108 


1.0732 x 10 6 


2 


Mean 


2.75 x 10~ 21 


0.91 


-0.00804 


5.31 


0.00108 


1.0733 x 10 6 




Std. Dev. 


±0.63 x 10" 21 


±0.28 


±0.52 


±0.46 


±0.000048 


±0.00019 x 10 6 




True 


8.54 x 10" 22 


4.66 


-0.800 


0.217 


6.15 


6.0000 x 10 b 




Recovered 


1.14 x 10" 21 


0.927 


-0.562 


5.45 


0.5* 


5.9991 x 10 5 


3 


Mean 


1.16 x 10~ 21 


0.962 


-0.488 


5.48 


0.24 


5.9992 x 10 5 




Std. Dev. 


±0.18 x 10~ 21 


±0.16 


±0.16 


±0.18 


±0.14 


±0.00076 x 10 5 



Table 3. As Table[T] but now for the three signals in the challenge data set. 



between the true and recovered parameters is broadly consistent with our expectations from 
the analysis of the training data, and once again the true parameters lie within a few standard 
deviations of the means in most case^jj] One surprise is that for source 1, the true value of 
/max is considerably below Nyquist, while our analysis recovered a flat distribution indicative 
of a value above Nyquist. We have carried out further checks, but there is no indication that 
the data favours a particular value of / max in this case. This is similar to what we observed 
for source 4 in the training data, which had a similar value of / max . However, in that case 
there was a narrow peak in the posterior in the vicinity of the true value, but no such peak was 
evident here. 

In Figure|4]we show the true waveforms, in the A channel, for the three sources present 
in the data set, and overlay them with the waveforms corresponding to each of the modes 
that we included in our submission. We see clearly that, in all three cases, the recovered 
waveforms reproduce the true waveforms well, with small differences that arise from noise 
fluctuations in the detector. As for the training data, we expect the recovered SNR to be a 
better measure of the quality of the solution than the values of the recovered parameters. In 
Table[4]we list the true SNRs in the A and E channels, and the SNRs recovered by each of the 
solutions included in the submission. As we found in the case of the training data, we have 
clearly recovered all of the SNR of the true signal, and the residual after subtraction would be 
at the level of instrumental noise. We thus regard the search as a success. 

In Figure [5] we show ID and 2D marginalised posterior probability distributions for the 
second of the sources present in the challenge data set. Comparing to the posteriors for the 
training source shown in Figure [3] we see that, in this case, / max can be measured very well, 
since it is far below Nyquist, but similar degeneracies in sky position and time of coalesence 
are clearly evident. In this case the sky position degeneracies are particularly severe since the 
source is at very low frequency. 

3.2. Model selection using Bayesian evidence 

The MLDC search demonstrated that MultiNest is able to correctly identify and 
characterize bursts from cosmic string cusps in LISA data. However, cosmic string cusps 

+ For source 1 , t c is quite far from the true value by this measure. However, this is due to the correlated degeneracies 
between t c , /3g and 4>s- The true time of coalesence was consistent with the full recovered posterior distribution. 
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Figure 4. Plots comparing the true waveforms with the waveforms for the various modes 
recovered by MULTlNEST and included in our submission for the analysis of the challenge 
data set. 
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Recovered SNR 


Source 


*c(s) 


Channel 


True SNR 


Mode 1 


Mode 2 


Mode 3 


1 


1.6 x 10 6 


A 
E 


41.0 
14.5 


41.2 
14.5 


41.0 
14.6 


N/A 


2 


1.1 x 10 6 


A 
E 


30.7 
13.9 


30.7 
13.9 


N/A 


N/A 




6 x 10 5 


A 


18.8 


18.9 


18.5 


18.4 


3 


E 


36.9 


36.7 


37.1 


36.8 



Table 4. Signal-to-noise ratios for the three sources in the challenge data set. We show the true 
SNRs in the A and E channels, plus the SNRs for all of the modes we submitted for analysis. 




Figure 5. As Figure|3] but now for the second of the sources identified in the blind challenge 
data set. 
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produce a very particular waveform signature. LISA might also detect bursts of gravitational 
radiation from other sources, such as supernovae of hypermassive objects, or even unmodelled 
sources. A search of the LISA data using a cosmic string cusp model could conceivably detect 
these other types of burst. In order to make scientific statements it is important to be able 
to say what the most likely origin for an observed burst could be. The Bayesian evidence 
provides a tool for carrying out model selection of this type, although it relies on having two 
alternative hypotheses to compare. As an alternative to a cosmic string model, we consider 
a sine-Gaussian burst, which is a generic burst waveform commonly used in burst analysis 
for ground-based gravitational wave detectors ifTTl . We can use the evidence ratio to compare 
these two models when analysing a data stream containing a cosmic string burst or containing 
a sine-Gaussian burst. There are two stages to this analysis — 1) determining the signal-to- 
noise ratio at which we begin to be able to detect a signal of each type in a LISA data set; 2) 
determining the SNR at which the Bayesian evidence begins to favour the true model over the 
alternative. 



3.2.1. Detection threshold To determine the SNR needed for detection of a burst source 
with the MultiNest algorithm, we generated a sequence of data sets containing a particular 
source and Gaussian instrumental noise generated using the theoretical noise spectral densities 



given in Eq. ( 15 I, but varying the burst amplitude between the data sets to give a range of 



SNRs between 5 and 15. Note that here and in the following we quote total SNRs for the A 
and E channels combined. Each data set was searched with MultiNest and the global log- 
evidence of the data set, as computed by MultiNest, was recorded. This was repeated for 
8 different cosmic string sources and 9 different sine-Gaussian sources. The parameters for 
the cosmic string sources were taken from the parameters of the five training and three blind 
sources injected in the MLDC data sets, as these covered the parameter space of possible 
signals nicely. The important parameters of the sine-Gaussian model are the frequency and 
the burst width, so we considered three different values of each, and constructed waveforms 
for the nine possible combinations, while choosing a random value for the sky position and 
waveform polarisation in each case. 

In Figure|6]we show the global log-evidence of the data as a function of the SNR of the 
injected signal for the cosmic string burst sources. Figure [7] shows the corresponding results 
for the sine-Gaussian bursts. The detection threshold for cosmic string bursts appears to be at 
an SNR of 7-11 depending on the source parameters. Most of the sources have a significant 
log-evidence (greater than 3) at an SNR of 7, but the SNR required for detection is higher 
for the first training source and the second blind source. These are the two sources with the 
lowest values of the break frequency / max , which makes the waveforms much smoother and 
simpler in the time domain, as indicated in Fig. [2] This may explain why it is more difficult to 
distinguish them from noise. The MLDC data sets had a prior for the source SNR, in a single 
Michelson channel, of p € [10, 100], so these results suggest we would be able to detect any 
source drawn from the MLDC prior. 

The sine-Gaussian bursts require a slightly higher SNR for detection, of 8-11, but this 
increase in threshold is only of the order of 1 in SNR. The sine-Gaussian signals are in general 
much simpler in form than the cosmic strings and therefore it is perhaps unsurprising that it 
requires a higher SNR to distinguish them from noise. In this case, it was the sources with 
highest frequency, f c = 0.49Hz, that were most difficult to detect. However, since the Nyquist 
frequency of the data sets was 0.5Hz, the difficulty of detection may have arisen because the 
signal was partially out of band. 
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Bayesian Evidence Ratio as a Function of SNR for a Cosmic String Signal 
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Figure 6. The Bayesian log-evidence for the data set as a function of the SNR for eight 
different cosmic string cusp signals. These signals had the same parameters, other than 
amplitude, as the various sources used in the training and blind data sets of the MLDC. The 
labels in the key are consistent with labels in Tables[T]andr3] 



Bayesian Evidence Ratio as a Function of SNR for a Sine Gaussian Signal 
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Figure 7. As Fig. [6] but now showing results for nine different sine-Gaussian signals, with 
frequency, /, and width, Q, as shown. 
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Figure 8. Typical example of confusion when searching a cosmic string data set with the 
wrong model. The plot shows a comparison of the injected cosmic string signal to the best- 
fit signals found by MULTlNEST using the cosmic string model as templates and using the 
sine-Gaussian model as templates. 



3.2.2. Model selection To explore model selection, we used MultiNest to search the 
same data sets described above, but now using the alternative model, i.e., we searched the 
cosmic string data sets using the sine-Gaussian likelihood and vice versa. For sufficiently 
high SNR, in both cases the alternative model was able to successfully detect a signal in the 
data set. Typically, the best-fit sine-Gaussian signal to a cosmic string source has low Q and 
a frequency that matches the two lobes of the cosmic string burst. The parameter Q sets the 
number of cycles in the sine-Gaussian wave packet, and so a sine-Gaussain with Q ~ 2 most 
closely resembles a cosmic string event, which typically has two cycles. This is illustrated for 
a typical case in Fig. [8] Similarly, the best-fit cosmic string source to a sine-Gaussian signal 
matches the central two peaks of the sine-Gaussian waveform as well as possible. A typical 
case is shown in Fig. [9] 

From the results of these searches it is possible to construct the evidence ratio of the two 
models for each of the data sets. In Fig.[l0]we show the ratio of the Bayesian evidence for 
the cosmic string model to that of the sine-Gaussian model when searching the cosmic string 
data sets. We see that the evidence ratio starts to significantly favour the true model, i.e., the 
cosmic string, at an injected SNR of ~ 7, which is the point at which we first start to be able 
to detect the cosmic string burst at all. For the two low frequency sources, training source 
1 and blind source 2, the evidence ratio only starts to favour the true model at SNR^ 11, 
but again this is the point at which we are first able to detect the source. We conclude that 
when a cosmic string burst is loud enough to be detected, then the evidence clearly favours 
the interpretation of the event as a cosmic string burst. 



In Fig. 1 1 we show the ratio of the evidence of the sine-Gaussian model to that of the 



cosmic string model when searching the data sets containing sine-Gaussian signals. The 
conclusions are broadly the same as for the cosmic string case. We require a slightly higher 
SNR in order to correctly choose the sine-Gaussian model, but this just reflects the fact that 
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Comparison of Recovered Signals for a Sine Gaussian 
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Figure 9. As Fig. [8] but we now show a typical confusion example for searches of the 
sine-Gaussian data sets. We compare the injected sine-Gaussian signal to the best-fit signals 
recovered by MULTlNEST using both the cosmic string and the sine Gaussian models. 
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Bayesian Evidence Ratio as a Function of SNR for a Cosmic String Signal 
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Figure 10. The ratio of the Bayesian evidence for the cosmic string model to that of the sine- 
Gaussian model when searching a data set containing a cosmic string burst source. We show 
the Bayesian evidence ratio as a function of signal SNR for a variety of different cosmic string 
sources. 
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Bayesian Evidence Ratio as a Function of SNR for a Sine Gaussian Signal 
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Figure 11. As Fig. |10| but now showing the ratio of the Bayesian evidence for the sine- 
Gaussian model to that of the cosmic string model when searching a data set containing a 
sine-Gaussian burst source. 



we need a somewhat higher SNR to detect the sine-Gaussians in the first place. The only case 
for which the evidence ratio does not begin to favour the sine-Gaussian model at the point 
where the source becomes detectable is the case with / = 0.002Hz and Q — 2. This is a 
sine-Gaussian signal with only two smooth oscillations, and so it does look rather like a low 
frequency cosmic string event. Even in that case, the evidence begins clearly to favour the 
correct model for SNRs of ~ 13 and higher. 

4. Discussion 

We have considered the use of the multi-modal nested sampling algorithm MultiNest for 
detection and characterisation of cosmic string burst sources in LISA data. As a search 
tool, the algorithm was able successfully to find the three cosmic string bursts that were 
present in the MLDC challenge data set. These sources, and the five sources in the MLDC 
training data, were correctly identified in the sense that the full signal-to-noise ratio of the 
injected source was recovered, and a posterior distribution for the parameters obtained. The 
maximum likelihood and maximum a-posteriori parameters were not particularly close to the 
true parameters of the injected signals, but this was a consequence of the intrinsic degeneracies 
in the cosmic string model parameter space and in all cases the true parameters were consistent 
with the recovered posterior distributions. 

In controlled studies, we found that the SNR threshold required for detection of the 
cosmic string bursts was ~ 7-11, depending on the burst parameters. Bursts with a low 
break- frequency require a higher SNR to detect than those with high break frequencies. We 
also explored the detection of sine-Gaussian bursts and in that case the SNR required for 
detection was slightly higher, being typically 8-11, with sources having frequency close to 
Nyquist being more difficult to detect. 

MultiNest is designed to evaluate the evidence of the data under a certain hypothesis, 
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and this can be used to compare possible models for the burst sources. LISA may detect 
bursts from several different sources, and it is important for scientific interpretation that the 
nature of the burst be correctly identified. We used the Bayesian evidence as a tool to choose 
between two different models for a LISA burst source — the cosmic string model and the 
sine-Gaussian model, which was chosen to represent a generic burst. The Bayesian evidence 
works very well as a discriminator between these two models. The evidence ratio begins to 
clearly favour the correct model over the alternative at the same SNR that the sources become 
loud enough to detect in the first place. 

The usefulness of MultiNest as a search tool in this problem is a further illustration 
of the potential utility of this algorithm for LISA data analysis, as previously demonstrated 
in a search for non-spinning SMBH binaries |7|. Other algorithms based on Markov Chain 
Monte Carlo techniques have also been applied to the search for cosmic strings [24]. Both 
approaches performed equally well as search tools in the last round of the MLDC. We are 
now exploring the application of MultiNest to searches for other LISA sources, including 
spinning massive black hole binaries and extreme-mass-ratio inspirals. MultiNest was not 
designed primarily as a search algorithm, but as a tool for evidence evaluation and this work 
has demonstrated the utility of the Bayesian evidence as a tool for model selection in a LISA 
context. Other problems where the evidence ratio approach could be applied include choosing 
between relativity and alternative theories of gravity as explanations for the gravitational 
waves observed by LISA, or choosing between different models for a gravitational wave 
background present in the LISA data set. The Bayesian evidence was previously used in a 
LIGO context as a tool to choose between alternative theories of gravity [15] and in a LISA 
context to distinguish a data set containing a source from one containing purely instrumental 
noise [16|. MultiNest provides a more efficient way to compute the evidence and this 
should be explored in more detail in the future. 

In the context of interpretation of LISA burst events, what we have considered here is 
only part of the picture. We have shown that we are able to correctly choose between two 
particular models for a burst, and this can easily be extended to include other burst models. 
However, LISA might also detect bursts from unmodelled sources. In that case, algorithms 
such as MultiNest which rely on matched filtering would find the best fit parameters within 
the model space, but a higher intrinsic SNR of the source would be required for detection. 
In such a situation, we would like to be able to say that the source was probably not from a 
model of particular type, e.g., not a cosmic string burst. There are several clues which would 
provide an indication that this was the case. The sine-Gaussian model is sufficiently generic 
that we would expect it, in general, to provide a better match to unmodelled bursts than the 
cosmic string model, which has a very specific form. Therefore, we could say that if the 
evidence ratio favoured the cosmic string model over the sine-Gaussian model it was highly 
likely that the burst was in fact a cosmic string and not something else. Similarly, if we found 
that several of the alternative models had almost equal evidence, but the SNR was quite high, 
it would be indicative that the burst was not described by any of the models. We have seen that 
at relatively moderate SNRs, when the signal is described by one of the models, the evidence 
clearly favours the true model over an alternative. If we found that two models gave almost 
equally good descriptions of the source, it would suggest that the burst was not fully described 
by either of them. A third clue would come from the shape of the posterior for the source 
parameters. The cosmic string waveform space contains many degeneracies, but these can be 
characterised theoretically for a given choice of source parameters. If the signal was not from 
a cosmic string, we might find that the structure of the posterior was modified. Finally, some 
techniques have been developed for the Bayesian reconstruction of generic bursts [29 1 which 
could also be applied in a LISA context. The usefulness of these various approaches can 
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be explored further by analysing data sets into which numerical supernovae burst waveforms 
have been injected. While the necessary mass of the progenitor is probably unphysically 
high for a supernova to produce a burst in the LISA frequency band, such waveforms provide 
examples of unmodelled burst signals on which to test analysis techniques. The final LISA 
analysis will employ a family of burst models to characterize any detected events. The work 
described here demonstrates that the Bayesian evidence will be a useful tool for choosing 
between such models, and MultiNest is a useful tool for computing those evidences. 
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